C  The following routine calculates and stores logs of factorials.
      SUBROUTINE factor(fl, nfl)
      integer nfl, i
      real*8 fl(0:nfl)
C
C ****  COMPUTE AND STORE, FL(I) = LOG((I)!)
C
      fl(0)=0d0
      do i = 1, nfl
         fl(i) = fl(i-1) + log(dfloat(i))
      end do 
      RETURN
      END
c------------------------------------------------------------------
      FUNCTION CPLM (L, M, X)
C     
C     CALCULATE ASSOCIATED LEGENDRE POLY FROM A RECURSION RELATION
c     for a complex argument
C     IN MESSIAH. NOTE: M CAN'T BE NEGATIVE
C     by Dima Konovalov dec.92 using the original code PLM

      complex*16 one, two, CPLM, x, A, B, DM, DF, OMX, DLL
      IMPLICIT integer (i-n)
      one = (1d0, 0d0)
      two = (2d0, 0d0)
c$$$      IF(ABS(X) .GT. 1.000001D0) THEN
c$$$         print*, ' WARNING : ERROR IN CPLM   abs(X) MUST BE less than 1',
c$$$     >      '  HERE X=', x
c$$$         stop ' Warning 1'
c$$$      END IF
      IF(L.LT.0) WRITE(6, 1) L
 1    FORMAT(' ? L IN CPLM = ', I6/)
      IF(M.LT.0) WRITE(6, 2) M
 2    FORMAT(' ? M IN CPLM = ', I6/)
      IF(M.GT.L) WRITE(6, 3) L, M
 3    FORMAT(' ? L = ', I6, ' & M = ', I6, ' IN CPLM'/)
C*********
      IF(L.NE.0) GO TO 4
      CPLM = one
      RETURN
C*********
4     IF(L.NE.1) GO TO 6
      IF(M.NE.0) GO TO 5
      CPLM = X
      RETURN
C*********
5     CPLM = SQRT(one-X * X)
      RETURN
C*********
C HERE  L.GE.2
C
6     IF(M.NE.0) GO TO 7
      A = one
      B = X
      DM = cmplx(M, 0)
      GO TO 10
C*********
C    CALCULATES CPLM(M, M, X)
C
7     DF = one
      ITMONE = 2 * M-1
      DO LL = 1, ITMONE, 2
         DLL = CMPLX(LL, 0)
         DF = DF * DLL
      end do 
      OMX = SQRT(one-X * X)**M
      A = OMX * DF
      IF(L.NE.M) GO TO 9
      CPLM = A
      RETURN
C**********
9     DM = cmplx(M, 0)
      B = A * (two * DM+one) * X
      IF(L.NE.M+1) GO TO 10
      CPLM = B
      RETURN
C**********
 10   MPTWO = M+2
      DO LL = MPTWO, L
         DLL = cmplx(LL, 0)
         CPLM = ((two * DLL-one) * X * B-(DLL+DM-one) * A)/(DLL-DM)
         A = B
         B = CPLM
      end do 
      RETURN
      END
C
C***
C
      FUNCTION PLM(L, M, X)
C
C     CALCULATE ASSOCIATED LEGENDRE POLY FROM A RECURSION RELATION
C     IN MESSIAH. NOTE: M CAN'T BE NEGATIVE
C

      IMPLICIT REAL*8 (A-H, O-Z)
      IMPLICIT integer (i-n)
      IF(ABS(X).GT.1.0D0) THEN
                             WRITE(6, 256) X
256   FORMAT(' WARNING : ERROR IN PLM   X MUST BE BETWEEN [-1, 1]',
     *        '   HERE X=', G22.16)
                            CALL EXIT
                        END IF
      IF(L.LT.0) WRITE(6, 1) L
1     FORMAT(' ? L IN PLM = ', I6/)
      IF(M.LT.0) WRITE(6, 2) M
2     FORMAT(' ? M IN PLM = ', I6/)
      IF(M.GT.L) WRITE(6, 3) L, M
3     FORMAT(' ? L = ', I6, ' & M = ', I6, ' IN PLM'/)
C*********
      IF(L.NE.0) GO TO 4
      PLM=1.0D0
      RETURN
C*********
4     IF(L.NE.1) GO TO 6
      IF(M.NE.0) GO TO 5
      PLM=X
      RETURN
C*********
5     PLM=DSQRT(1.0D0-X*X)
      RETURN
C*********
C HERE  L.GE.2
C
6     IF(M.NE.0) GO TO 7
      A=1.0D0
      B=X
      DM=M
      GO TO 10
C*********
C    CALCULATES PLM(M, M, X)
C
7     DF=1.0D0
      ITMONE=2*M-1
      DO 8 LL=1, ITMONE, 2
      DLL=LL
8     DF=DF*DLL
      OMX=DSQRT(1.0D0-X*X)**M
      A=OMX*DF
      IF(L.NE.M) GO TO 9
      PLM=A
      RETURN
C**********
9     DM=M
      B=A*(2.0D0*DM+1.D0)*X
      IF(L.NE.M+1) GO TO 10
      PLM=B
      RETURN
C**********
10      MPTWO=M+2
       DO 11 LL=MPTWO, L
      DLL = DBLE(LL)
      PLM=((2.0D0*DLL-1.0D0)*X*B-(DLL+DM-1.0D0)*A)/(DLL-DM)
      A=B
      B=PLM
11    CONTINUE
      RETURN
      END

